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Abstract: We present a hierarchical approach for enhancing the robustness 
of numerical solvers for modelling radiative MHD flows in multi-dimensions. 

This approach is based on clustering the entries of the global Jacobian in 
a hierarchical manner that enables employing a variety of solution procedures 
ranging from a purely explicit time-stepping up to fully implicit schemes. 
A gradual coupling of the radiative MHD equation with the radiative transfer 
equation in higher dimensions is possible. 

Using this approach, it is possible to follow the evolution of strongly time- 
dependent flows with low/high accuracies and with efficiency comparable to 
explicit methods, as well as searching quasi-stationary solutions for highly 
viscous flows. 

In particular, it is shown that the hierarchical approach is capable of modelling 
the formation of jets in active galactic nuclei and reproduce the corresponding 
spectral energy distribution with a reasonable accuracy. 

Key words: Methods: numerical - hydrodynamics - MHD - radiative 
transfer 



1 Introduction 



Within the last two decades, a tremendous progress has been made in both 
computational fluid dynamics (CFD) algorithms and the computer hardware 
technologies. The computing speed and memory capacity of computers have 
increased exponentially during this period. Similarly is in astrophysical fluid 
dynamics (AFD), which is a rapidly growing research field, and in which mod- 
ern numerical methods are extensively used to model the evolution of rather 
complicated flows. Unlike CFD, in which implicit methods are frequently used, 
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the majority of the methods u sed in AFP are expli cit. Se veral of t hem be - 
came very po pular, e.g., ZEUS jST92ft . NIRVANA+ ilZIQSl) . FLASH iFROft) . 
VAC l|T098|) . THARM l|GA03l) . The popularity of explicit methods arises 
from their being easy to construct, vectorizable, parallelizable and even more 
efficient as long as dynamical evolutions of compressible flows are concerned. 
Specifically, for modeling the dynamical evolution of HD-fiows in two and 
three dimensions explicit methods are hig hly superior to -date. For modelling 
relativistic flows, Koide and collaborators (|K099t iKOnl and l|ME0ltlKO99l 
have developed pioneering general relativistic MHD solvers. A rather com- 
ple te review of n umerical approaches for relativistic fluid dynamics is given 
in l|MA99tlFO00|) . A ZEUS-like sch eme fo r general relativistic MHD has also 
been developed and is described in 

These methods, however, are numerically stable as far as the Courant- 
Friedrich-Levy number is smaller than unity. The corresponding time step 
size decreases dramatically with the incorporation of real astrophysical effects. 
Specifically, they may even stagnate if self-gravity, radiative and chemical ef- 
fects are included. Moreover, explicit methods break down if the flow is weakly 
or strongly incompressible, and if the domain of calculations is subdivided into 
a strongly stretched mesh. In an attempt to enhance their robustness, several 
alternatives have been sug gested, such a s semi-explicit, semi- implicit or even 
implicit-explicit methods l|KL89t |T098). Nevertheless, their rather limited 
range of applications has lead to the fact that most of the interesting astro- 
physical problems remained, indeed, not really solved. A simple example is 
the evolution of a steady turbulent accretion disk. It was found by Balbus 
& Hawley (1991) that weak magnetic fields in accretion disks are amplified, 
generate turbulence, which in turn redistribute the angular momentum in the 
disk. However, whether this instability leads to the long-sought global steady 
accretion rate, or is it just a transient phenomenon in which the generation of 
turbulence is subsequently suppressed by dynamo action are not at all clear. 
Other notable phenomena are the formation and acceleration of the observed 
superluminal jets in quasars and in microquasars, the origin of the quasi- 
periodic oscillation in low mass X-ray binaries or the progenitors of gamma 
ray burst are still spectacular. 

Explicit methods rely on time-extrapolation procedures for advancing the 
solution in time. However, in order to provide physically consistent solutions, it 
is necessary that these procedures are numerically stable. The usual approach 
for examining the stability of nu merical methods is to perform the so called 
von Neumann analysis (sec HI90, for further details). This yields the so called 
Courant-Friedrich-Levy condition (CFL) which is known to limit the range of 
application and severely affects the robustness of explicit methods. In particu- 
lar, equations corresponding to physical processes occurring on much shorter 
time scales than the hydro-time scale (e.g., radiation, self-gravitation and 
chemical reactions) cannot be followed explicitly. Furthermore, these meth- 
ods are not suited for searching solutions that correspond to evolutionary 
phases occurring on time scales much longer than the hydro-time scale. Us- 
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ing high performance computers to perform a large number of explicit time 
steps may lead to accumulation of round-off errors that can easily distort the 
propagation of information from the boundaries and cause divergence of the 
solution procedure, especially if Neumann type conditions are imposed at the 
boundaries. In contrast to explicit methods, implicit methods are based on 
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Fig. 1. A schematic description of the time step size and the computational costs 
versus the band width M of the Jacobian. N is the number of unknowns. Explicit 
methods correspond to M = 1 and large 1/St. They require minimum computational 
costs (CC). Large time steps (i.e., small 1/St) can be achieved using strongly implicit 
methods. These methods generally rely on the inversion of matrices with large band 
width, hence computationally expensive, and, in most cases, are inefficient. 

solving a matrix equation of the form Ax = 6, where A is the Jacobian matrix 
corresponding to the system of equations to be solved, b is the right hand side 
vector of known quantities, and x is the solution vector sought. These methods 
have two major drawbacks. First, constructing the matrix A is difficult, time 
consuming, and may considerably influence the robustness of the method. 
Second, the inversion procedure must be stable and efficient. In general, con- 
servative discretization of the MHD equations give rise to sparse matrices, or 
even to narrow band matrices. Therefore, any efficient matrix inversion pro- 
cedure must take the advantage of A being sparse. Inverting A directly by 
using Gaussian elimination requires N 3 algebraic operations, where N is the 
number of unknowns. If the flow is multi-dimensional and a high spatial reso- 
lution is required, the number of operations can be prohibitive even on modern 
supercomputers. Krylov Sub- Iterative Methods (KSIMs), on the other hand, 
are most suited for sparse matrices and avoid the fill-in procedure. In the 
latter case, A is not directly involved in the process, but rather its multiplica- 
tion with a vector. The convergence rate of KSIMs has been found to depend 
strongly on the proper choice of the pre-conditioncr. For advection-dominated 
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flows, incomplete factorization such as ILU, IC and LQ, approximate factor- 
ization, ADI, line G auss-Se idel are only a small sub-set of possible sequential 
pre-conditioners fsee lSAOol and the references therein). Another powerful way 
of accelerating relaxation techniqu es is to use th e multi-grid method as a di- 
rect solver or as a pre- conditioner <|BR01HTR0i|) . For parallel computations, 
Red-Black ordering in combination with GRMES and Bi-CGSTAB as wel l as 
domain decomposition are among the popular pre-conditioners fsee lDQ98l for 
further discussion.) 

Towards studying the jet-disk-BH connection in AGNs and 11— QSOs a se- 
ries of multi-dimensional cal culations have been performed (e.g.. lOU97tlUC99T : 
\MFM iHTTjj lHIT.1,4 iHTJJfih . Specifically, these studies revealed that: 

1. Counter-rotating disks with respect to the BH-spin generate jets that 
propagate approximately twice as fast as in the co-rotating case. 

2. Jets formed are found to be relatively slow, i.e., the corresponding- 
ly—factors did not reach the desired large values. This was found in both 
cases: when the spins of the disk and the BH are parallel and when they 
are anti-parallel. Moreover, disks surrounding Kerr BHs have been ver- 
ified to produce jets that are more powerful than in the Schwarzschild 
case. These jets are driven primarily by strong MFs that are created by 
the frame dragging effect. 

3. Large factors are obtainable if the Alfven speed due to the PMF is 
equal to or even larger than the local escape velocity (see ME03, and the 
references therein). 

4. Poliodal magnetic fields may extract rotational energy from the disk 
plasma, and from a geometrically thin super-Keplerian layer between the 
disk and the overlying corona. The outflowing plasma in this layer is dis- 
sipative, two-temperature, virial-hot, advective and electron-proton dom- 
inated. The innermost part of the disk in this model is turbulent-free, 
sub-Keplerian rotating and advective-dominated. This part ceases to ra- 
diate as a standard disk, and most of the accretion energy is converted 
into magnetic and kinetic energies that go into powering the jet. 

Nevertheless, jet-structures, their formation, acceleration, their linkage to 
the accretion phenomena and the nature of their plasma are still a matter 
of debate. Furthermore, the flood of observational data makes it even more 
essential than ever to perform sophisticated numerical calculations to gain a 
more precise insight of their evolution. 

In this paper we focus on the architecture of the global solution procedure 
rather than on local details, such as order of accuracies, physical consistency, 
types of advection schemes or fulfilling the solenoidal condition. Specifically, 
we discuss strategies for enhancing the robustness of solvers through con- 
structing various pre-conditionings to implement a variety of solution methods 
in arbitrary dimensions. Special attention is given to radiative MHD solvers 
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and their possible coupling with the radiative transfer equation in higher di- 
mensions. 



2 The governing equations 

2.1 The 3D axi-symmetric radiative MHD equations 

Spherical geometry is the most appropriate geometry for capturing flow con- 
figurations in the vicinity of black holes. Taking into account the perfect axi- 
symmctry of black holes, and that their gravitational pull dominates the forces 
exerting on the surrounding flows, we conclude that axi-symmetry is a rea- 
sonable assumption that may characterize accretion flows in their vicintics. 
Moreover, in applying spherical geometry the transformation 8 = tt/2 — 9 has 
been used 1 . We note that the dynamical time scale near the event horizon is 
extremely short, therefor giving rise to multi-component flows, such as elec- 
tron and ion plasmas. 

In the following we describe the set of radiative MHD equations, and list the 
scaling variable that may be used for transforming them into non-dimensional 
form (see Table 1). 

• Continuity equation: 

g+V-pV-0 (1) 

• Radial momentum equation: 

dm „ dP (V? + V v 2 ) d$ , dE ^ r ^ r . , 

- d J + V - mV = ^+P r +^ + AFLD ^r + F £ + Qvis (2) 

• Vertical momentum equation: 

dn „ Tr dP lr9 „ d<P dE „ 

at + ' = W ~ pv * ' + p de + Afld ~ae + + Qvis (3) 



+ V-£V = FZ + Q* is (4) 



• Angular momentum equation: 

at 

• Internal equation of the ions: 

dE d 

+ V • EfV = -( 7 - l)£ d V -V + $- yli-e + V • < ond VTi (5) 

• Internal equation of the electrons: 

df d 

-^ r +\7-S*V = -( 7 -l)f c d V-y+yl i _ c -/lB-/lc-^Sy„+V-C nd VT c (6) 



1 This transformation allows simple analogy with and into cylindrical coordinates. 
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• Equation of the zero moment of the radiation field: 

f) E 

— + V • EV + V • [A FLD V£] - A B + A c + A Syn (7) 

• The induction equation: 

dB 

— = V x (V x B + a dyn B - */ mag V x B). (8) 

• Gravitational potential: the Poisson equation: 

Atp = inGp, (9) 
where ip is the gravitational potential and G is the gravitational constant. 
In Table (2) we list part of the variables used and their definitions. Further, 



Scaling variables: 

Mass: M = 3 x 1O 8 M 

Accretion rate: M = 10 _1 Msdd 

Distance: R = R in = 3R S , where R s = 2GM /c 2 

Temperature: f = 5 x W 7 K 

Velocities: V = V s = [7^ ga st//ii] 1/2 , m = 1.23 

Ang. Velocity: % = Vk cp = {GM/R) 1/2 
Magnetic Fields: B = Vs/V^jFp 

Density: p = M/{H d R~ ^Vs) = 2.5 x 10~ 12 gem' 3 

Table 1. Scaling variables for reformulating the MHD equations in non-dimensions. 



the subscripts "i" and "e" correspond to ion and electron plasmas, where 
7 = 5/3, jUj = 1.23 and /i G = 1.14 are used, oyyn, Wg correspond to the 
a— dynamo and the magnetic diffusivity, respectively. The radiative diffusion 
coefficient Afld is a radiative flux limiter which forces the radiative flux to 
adopt the correct form in optically thin and thick regions, i.e., 

and provides a smooth matching in the transition regions. Here \ — p(K a bs+o") 
and n = VE/|VE|, where K a bs and a are the absorption and scattering coef- 
ficients. Ab, Ai-e, Ac, vl S yn correspond to Bremsstrahlung cooling, Coulomb 
coupling bet ween th e ions and electrons, Compton and synchrotron coolings, 
respectively l|RY79l) . These processes read: 

(T- — T 1 

yli_ c = 5.94 x 10- 3 nin e ck y 1 - . e IM 

e 

A B = AacK &hsP {T 4 - E)/N, (11) 
A c - 4 CT n e c(-^)(T e - T rad )E/M, 
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where N — [(7 — 1) / -f](V 2 V v / R) is a normalization quantity. n e , n\ are the 
electron- and ion- number densities. E is the density of the radiative energy, i.e., 
the zero-moment of the radiative field. The radiative temperature is defined 
as T ra d = E 1 / 4 . The Lorenz forces acting on charged plasma in the MHD 
approximation read: 



Fl = B p -VB = B t — (rcos6B T ) + -J- — (r cos 6 B T ). 



_ (rcos ^ T ) + __ 
The turbulent-diffusive terms read: 

T, 



Qvis ~ 


f 


!<-^>+ 


f 

r cos 6 


d 




f 


|<^)+ 


f 

rcos£ 


d 


Qt = 

^ VIS 


f 




f 

r cos( 


d 

we 



(cos er I( 



r 



(cos 0T M )+T^ tan 6 (13) 
(rcos 2 flT ev ), 



where 



ar 3 r or r cos fl oO 

Tee = 2r ? — a# + q (~ 2" o + 77757T cos^e 

r afl r 2> r z or r cos afl 

T 00 = 2rK- tanfl + ^ - + -J— |-(cos0^))) (14) 

r r 3 or r cos fl afl 

cosfl 9 
r de cosfl 

a V e ,,ldV t 
ar r r afl 

2.2 The isotropic radiation transfer equation: The Kompaneets 
equation 

Compton up-scattering of soft photons is most efficient in unsaturated Comp- 
tonization regions where the Compton-Y parameter is of order unity. This 
parameter acquires large values in optically thick media, and small values 
in the corona, implying that the corona-disk interaction region and/or the 
innermost region of the disk arc most appropriate for this process to op- 
erate efficiently. As a consequence, Comptonization in accretion flows is in- 
trinsically two-dimensional, and therefore requires a multi-dimensional treat- 
ment. So far, Comptonization has been considered under strong assumptions 
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Symbols: 


V 


= (Vr, V0, Vy) 


velocity field 


li 


— (ft u n n„\ — (ft 


magnetic licit! 


V 


- (A. IS-) 

\ fir ' r ae 1 


gradient in spherical coordinates 


V- 


= 4 s r 2 + 1 „ LcosO 

r z r r cos # of 


divergence in spherical coordinates 


rpe,i 




electron and ion temperatures 




= 1Z e asp(Ti/fii +T e /fx e )) 


electron and ion pressure 




= P e 'V(7-l), 


electron and ion internal energies 




= 7.8 x Te /2 , 3.2 x Tf /2 


electron and ion conductivities 


(m, n, Q 




momentum 






turbulent and magnetic diffusivities 




= ^HD + ^MHD 


HD and MHD turbulent dissipation (see MM) . 



Table 2. Variables used and their definitions 



that allow separation of variables and lead to the separation of the Kom- 
paneets operator from the radiative transfer equation. Here, the radiative 
intensity is assumed to be time-independent, isotropic and the plasma is 
isothermal. In this case, the generation and Comptonization of photons can 
be de s cribed by a second order d ifferential equation in the frequency space 
(|TL72l lFE72t \k A 74 ISH74 IHTJ.T.^ . 

Different accretion models display different spectra. Therefore, it is essen- 
tial to perform a diagnostic study to analyze their consistency with observa- 
tions. This however requires solving the 7D radiation transfer equation: 



1 dl 
~c~dt 



n • VI = K v p{S v -I)-apI + J CI dQdE + e™ d , (15) 



where I = I(t, r, 9, if, </>, v) is the radiative intensity which depends on time 
t, the spherical coordinates (r,9,tp), two ordinates <f>) that determine the 
direction of the photons on the unit sphere, and on the frequency v. k v and a 
are the absorption and scattering coefficients. S v is a source function. 7i nt = 
j CI dfldE describes the scattering of photons through electrons, and C is 
the scattering kernel. £™ od is the modified synchrotron emission. 

To make the problem tractable, the following approximations have been 
performed: 

• The radiation field is axi-symmetric and isotropic, i.e., d/dip = and 

J = t- f idn « /. 

• The source function is represented by the modified black body function, 
i.e., 

S v = i?r d = , (16) 

1 + Jl + 



where B v is the normal Planck function fsee lRY7l^) . 

The thermal energy of the electrons is far below its corresponding rest 
mass energy, i.e., e = kT •> <C 1, and hv . < , <C 1. 
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Using the last approximation, / int can be expanded up to sec ond order in e 
which reduces it to the so-called Kompaneets operator l)PA80|) : 

Zint JCu = j i- (4fcT - hv)I+ ^H. (17) 

m c c ov m c c cw 

In this case, the radiative transfer equation with respect to a rest frame of 
reference reads: 

1 F)F* A 

-hsr + v • VE ^ = - A -( v • V ) E » + v • [— VE ^ 

c at Xv 

+ K 1/ p(S v -E v )+JC u +e™ d , (18) 

where E v = —f-J = E(t, r, 9, v),Xv = p( K ^ + (T ), A„ is the flux limited diffusion 
coefficient |LE81), which forces the radiative flux to adopt the correct form 
in optically thin and thick regions, i.e., 

X u may provide a smooth matching between these two extreme regimes. The 
above two different behaviour of the operator can be combined as follows: 

V ■ KV V E U ^ V • r) T VE u , (20) 

where r ) r = (l-aO^f^ + agj, a = er RFLD , and i? FL D = V E / p{n v S v + oE u ) 
\\\VM\ \. 

£ mod - n ^g) corresponds to the modified synchrotron emission of pho- 
tons by relativistic electrons gyrating around magnetic field lines, which reads: 

where f(= e'^^ ) is a switch on/off operator which bridges optically thin 
and thick media to synchrotron radiation, and v c is a critical frequency (see 
below) . 

An appropriate approximation for e v in optically thin medium reads 
(|MA9fJl : 

e v = 2.73 x 1(T 5 P v j( y> B, 6) ergs cm~ 3 s" 1 Hz" 1 , (21) 
where K 2 is the Bessel function of the second kind and I = ^| (1 + -^j^ + 

053\ P -1.89C 1/3 

£1/2 J e 

Here C = 2.38 x 10" 7 ' (v / B6 2 ) and 6 e = kT c /m c c 2 . 

Below a certain critical frequency v c , the media become self- absorbing to 
synchrotron emission. In this case, e™ od « e BB = To find v c , we use 

the local non-linear Newton iteration procedure applied to the equation 



e v dV= / e^dS. (22) 



J3C, 

'S 

Having obtained v c , the switch on/off operator £ can then be constructed. 
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Fig. 2. A schematic description of the hierarchical solution method. A cluster 
of coefficients is computed in the first stage, and a matrix-generator is created that 
allows using various solution procedures ranging from purely explicit to fully implicit. 
Interchange between solution methods is possible, as modifying, adding or removing 
entries is directly maintainable. 
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3 Solution methods 

3.1 Solving the radiative MHD equations 

The set of equations in conservative form may be written in the following 
vector form: 

<9q 

— + L TII F + Lg_ggG = f , (23) 

at 

where F and G are fluxes of q, and L r rr , Lg t gg are first and second order 
transport operators that describe advection-diffusion of the vector variables 
q in r and 9 directions, f corresponds to the vector of source functions. 
Adopting a five star staggered grid discretization, it is easy to verify that at 
each grid point the Eq. (23) acquires the following block matrix equation: 

^ + S*S qj _ hk + D*6q. hk + S r Sq j+lik 

+ ^dq^ + D e Sq jM + S 9 Sq hk+1 = RHS? k , (24) 

where the subscripts "j" and "k" denote the grid-numbering in the r and 
9 directions, respectively, and RHS n = [f — £ r . rr F — Lg,ggG] n . Underlines 
(overlines) mark the sub-diagonal (super-diagonal) block matrices in the cor- 
responding directions, and D r ' 6 are the diagonal block matrices. 
To outline the directional dependence of the block matrices, we re-write Eq. 
24 in a more compact form: 

S ^j.k+i _ 

+£ r <fy_i,k +Anod<ty,k +S r &?j+i,k = RHSf <k (25) 
+S e Sq ik _ 1 , 

where £> mo d = Sqj k /5t + D x + D y . Eq. (25) gives rise to at least four different 
types of solution procedures: 

1. Classical explicit methods are very special cases in which the sub- and 
super-diagonal block matrices together with D x and D y are neglected. 
The only matrix to be retained here is (l/St) x (the identity matrix), i.e., 
the first term on the LHS of Eq. 24. This yields the vector equation (see 
M5/Fig. 2): 

\^5q- hk = RHSf >k . (26) 

2. Semi-explicit methods are obtained by preserving the diagonal entries, 
dj,k, of the block diagonal matrix D mo d (see M4/Fig. 2). This method has 
been verified to be numerically stable even when large Courant-Friedrich- 
Levy (CFL) numbers are used. In particular, this method is absolutely 
stable if the flow is viscous-dominated. 
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Fig. 3. A schematic description of the hierarchical algorithm for solving the radia- 
tive MHD equations. Stage I corresponds to the implicit operator splitting approach 
(IOS), which is most appropriate for following the early time-dependent phases of 
the flow. The solution obtained can then be used as initial condition for Stage II, 
where the hydro-equations are solved as a single coupled system, followed by the 
magneto component, which is again solved as a single coupled system. Here, high 
spatial and temporal accuracies in combination with the prolongation/restriction 
strategy may be used. Similarly, the solution obtained in this stage may be used as 
starting solutions for Stage III, where steady solutions for the fully coupled set of 
equations consisting of the zero moment of the radiation field and the MHD equa- 
tions are sought. In this stage, pre-conditioned Krylov sub-iterative methods are 
considered to be robust and efficient. The very last stage, Stage IV, corresponds to 
the case where solutions for the internal energy equations weakly coupled with the 
5D radiative transfer equation are sought. 
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3. Semi-implicit methods are recovered when neglecting the sub- and super- 
diagonal block matrices only, but retaining the block diagonal matrices 
(see M3/Fig. 2). In this case the matrix equation reads: 

D mod 5q iik = RHSfa. (27) 

We note that inverting Z? mo d is a straightforward procedure, which can 
be maintained analytically or numerically. 

4. A fully implicit solution procedure requires retaining all the block ma- 
trices on the LHS of Eq. 25. This yields a global matrix that is highly 
sparse ( Ml/Fig . 2). In this case, the "Approximate Factorization Method" 
(-AFM : lBW78|i and the "Line Gauss-Seidel Relaxation Method" (-LGS: 
MA85) are considered to be efficient preconditionings for the set of radia- 
tive MHD-equations. 

In the case that only stationary solutions are sought, convergence to steady 
state can be accelera ted by adopting the so called the "Residual Smoothing 
Method" (see IHUJ81 and the references therein). 

This method is based on associating a time step size with the local CFL- 
number at each grid point. While this strategy is efficient at providing quasi- 
stationary solutions within a reasonable number of iterations, it is incapable 
at providing physically meaningful time scales for features that possess quasi- 
stationary behaviour. Here we suggest to use the obtained quasi-stationary 
solutions as initial configuration and re-start the calculations using a uniform 
and physically relevant time steps. 

3.2 The 5D axi-symmetric RT equation: method of solution 

Let CE = be the equivalent operator form of Equation 18 in the continuous 
space fic- £■£ consists of several terms, each of which requires a careful and 
different representation in the finite discretization space flh which is defined 
as [ti,t 2 ,...,t N ] <g) [rx,r 2 ,...,rj] ® [9 1} 6 2 , ... , Ok] ® [vi,v2,—,vm]- M, [r], [6] 
and [v] correspond to time, radius (spherical), latitude, and to the frequency 
intervals, respectively. 

In most astrophysical problems, radiative effects occur on relatively short 
time scales compared to the hydro- or magneto-hydrodynamical ones, for 
which the use of unconditionally implicit numerical solvers is essential. This 
requires however that all terms of Eq. (18) should be evaluated on the new 
time-level. The discretization used should assure that the resulting Jacobian 
ArBv — dCE/dE is diagonally dominant. Therefore, the following procedures 
are employed. 

• The advection term W -VE V is discretized using a second order up- winding. 
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Fig. 4. The problem of free-fall of gas onto a Schwarzschild black hole. The evolution 
of the CFL-mimber and the residual versus number of iterations are shown, using 
different solution procedures. The solution methods are: normal explicit (top/left), 
semi-explicit (middle/left), semi-explicit in combination with the residual smoothing 
strategy(bottom/left), semi-explicit using moderate CFL-numbers (top/right), semi- 
explicit method in which the time step size is taken to be a function of the maximum 
residual (middle/right), and finally the fully implicit method (bottom/right). The 
different forms of the semi-explicit method used here are stable and converges to 
the stationary solution, though at remarkably different rates. 
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Fig. 5. Free-fall of gas onto a black hole surrounded by a static cold disk. Top: the 
density distribution (red: large values, blue: low values, green: intermediate values). 
Middle: the temperature distribution (red: large values, blue: low values, gray: in- 
termediate values). The curved shock front, where the temperature attains maxima 
is obvious. Bottom: the distribution of the velocity field is shown. 
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• The second order diffusion term V ■ [A„V1?„] is discretized using second 
order central-difference scheme on a staggered grid 

• K, v contains advection and diffusion terms in the frequency space. Here 
up-winding discretization in the frequency space is used. 

Combining the contributions of all terms of Eq. (18), we obtain at each grid 
point the following equation: 

CT pncw i ~o r 77incw i Q@ pncw 

3- ^j-l.kan T & ^j+l.k.m — ^j.k-l.m 

i ~q® pncw i qv pncw ■ ~Q ly pnew 

+ ID 1 ' + D e + D»)E?™ m = RHS, (28) 

where S T = dCE„/dE^ 1Xm ,S'' = dCE u /dE i+1Xm , D r + D e + W = dCE u /dE iX 
S e = d£E u /dE hk _ hm , S 9 = 8£E u /dE Sik+ltm , T = 8CE v /dE 1Xm+1 and 
= dCE„/dEj k m _i. The terms 5 r ,£> r , and S correspond to the sub- 
diagonal, diagonal and super-diagonal entries of the Jacobian A r g„ in the 
radial direction respectively. A similar description applies to the 9— and 
^—directions. 

Thus, solving the equation at all grid points, is equivalent to solve matrix 
equation: A r8u E ncw = E old , or simply, Aq = b. 

This matrix is highly sparse, and pre-conditionings such as the Alternat- 
ing Direction Implicit (ADI) and the Approximate Factorization Methods 
(AFM) are considered to be efficient. However, ADI is not appropriate for 
searching steady solution in three or more dimensions, as it is numerically 
unstable in high dimensions l|FL88f) . Alternatively, we have tried the AFM as 
a pre-conditioner. However, it turns out that the AFM converges slower than 
our favorite iterative metho d: 'Bl ack- White-Brown' line Gauss-Seidel method 
(henceforth BWB-LGS, see 11100. for further details). The latter method pre- 
serves the diagonal dominance of A, and hence converges faster than AFM. 
It should be noted that the line Gauss-Seidel method in its classical form is 
not appropriate for vector and parallel machines, mainly because the vector- 
length is proportional to the number of unknowns in one direction. A reason- 
able way to extend the vector-length is to solve for all unknowns located on 
even-numbered grid points, and subsequently on odd-numbered grid points. 
The resulting vector-length in this case is proportional to the number of un- 
knowns in the plane under consideration, and therefore enabling enhancement 
efficiency when using vector or parallel machines. 

More specifically, in each plane we perform two sweeps: in the first sweep 
we consider the unknowns in the r — 9 plane, i.e., we solve the system of 
equations: 

S T SE^J Xm + (D* + D e + D v )5Ef^ m + S'SEf^ = RHS, 
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where j = 1 — > J and k runs over odd- numbered rows. In the second sweep, 
we solve: 

S r 5E™Z Km + (D* + D e + D")5E?g m + tfSEPg^ 

= RHS + S e 5Ef™ lui + S 5£j°k+i iin , 

where j = 1 — ► J and k here runs over even-numbered rows. Therefore, 
we actually perform 6-inversion procedures per each time step. Here the 3- 
dimensional problem is replaced by three one-dimensional problems that are 
solved iteratively to recover the solution of the original problem. The method 
is relatively efficient, as the overall number of arithmetic operations scales 
linearly with the number of grid points (~ 6 x 9 x N). 



Fig. 6. The distribution of the Bernoulli number of accretion flows around a su- 
permassive black hole. The Bernoulli number characterizes the energies of the flow 
in different regions. Gravitationally bound flows have negative total energy, whereas 
flows of positive total energy are gravitationally unbound, and potentially should 
expand to infinity. In this figure, the decrease of the Bernoulli number from large to 
low positive values is represented by yellow, green and red colours, whereas the blue 
colour corresponds to negative values. Obviously, gravitationally unbound blobs are 
formed in the vicinity of the black hole, which thereafter collimate under the action 
of magnetic fields to form the observed highly collimated jets. 
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Fig. 7. The initial distribution of the density in an accretion disk (M = 0.1 A^Edd) 
around a Schwarzschild black hole overlied by coronal plasmas (solid lines). The 
dashed lines correspond to the magnetic field lines threading the disk and the corona. 
The distance is given in units of 2.75 Rsch, where i?sch is the Schwarzschild radius. 
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Fig. 8. The evolution of the residual in the maximum norm versus the number 
of iteration. The adopted density and temperature profiles correspond to standard 
accretion disks (see Fig. 7). 



4 Validation and preliminary tests 

4.1 Free-fall of plasma onto a Schwarzschild black hole 



A centrifugally-unsupported gas around a spinless black hole is gravitationally 
bound, and therefore should fall-freely onto the black hole, provided that no 
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Fig. 9. A VLA image shows the active central engine of the giant elliptical galaxy 
M87 (top), and a NRAO radio image of the jet apparently emanating from within 
100 gravitational radii. Solid lines correspond to calculated profiles and the asterisks 
to observational data. The profiles 01 to 06 show the spectral energy distribution 
calculated using different magnetic field strengths, or different truncation radii, or 
high/low corona temperatures. In particular, the profile 07 corresponds to a model 
in which the toroidal magnetic field is set to vanish artificially, whereas the poloidal 
magnetic field is set to be in equipartition with the thermal energy of the electrons. 
The profile 08 is similar to 07, except that the toroidal magnetic field is allowed to 
develop and reach values beyond equipartition with respect to the thermal energy 
of the electrons in the transition layer between the disk and the overlying corona. 
The above spectral energy distribution has been obtained by solving the radiative 
transfer equation in 5-dimensions, taking into account the Kompaneets operator 
for consistently modelling Comptonization. 400 non-linearly distributed frequency 
points have been used to cover the frequency-space, and 125 x 40 finite volume cells 
to cover the spatial domain of the calculation. 
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other external forces oppose gravity. In this case, the radial distributions of 
the density and velocity far from the event horizon obey the power laws: r~ 3 / 2 
and r -1 / 2 , respectively. 

This physical problem is relevant for testing the flexibility of the hierarchical 
scenario at adopting various solution methods, and to test their capability to 
capture steady, oscillation-free and advection-dominated flows, even when a 
strongly stretched mesh distribution is used. 

The equations to be solved in this problem are the continuity, the radial and 
horizontal momentum equations, and the internal energy equation. The flow 
is assumed to be inviscid and adiabatic (7 = 5/3). The equations have been 
solved using a first order accurate advection scheme both in space and time. 
In carrying out these calculations, the following conditions/inputs have been 
taken into account: 

• The central object is a one solar-mass and non-rotating black hole. 

• The outer boundary is 100 times larger than the the inner radius, i.e., 
^out = 100 x i?j n , where i?; n is taken to be the radius of the last stable 
orbit 2 i?LS ■ To first order in V/ c, the flow at this radius can be still treated 
as non-relativistic, though the error can be as large as 30%. 

• Along the outer boundary, the density and temperature of the gas assume 
uniform distributions, and flow across this boundary with the free-fall ve- 
locity. Symmetry boundary conditions along the equator, and asymmetry 
boundary conditions along the axis of rotation have been imposed. Along 
the inner boundary, we have imposed non-reflecting and outflow condi- 
tions. This means that up-stream conditions are imposed, which forbid 
information exterior to the boundary to penetrate into the domain of cal- 
culations. In particular, the actual values of the density, temperature and 
momentum in the ghost zone r are erased and replaced by the correspond- 
ing values in the last zone, i.e, the zone between R ln and R ln + AR. In 
the case that second order viscous operators are considered, care has been 
taken to assure that their first order derivatives across R m are vanished. 

The above set of equations are solved in the first quadrant [1 < r < 100] x [0 < 
9 < n/2], where 200 strongly stretched finite volume cells in the radial direc- 
tion and 60 in the horizontal direction are used. 

In Fig. 4, we show the evolutions of the CFL-number and the residual as 
function of the number of iteration which has been obtained using various 
numerical approaches. The results show that the convergence of the explicit 
and semi-explicit methods are rather slow when a relatively small time step 
size is used. This implies that the amplitude-limited oscillations are strongly 
time-dependent that may result from geometric compression. Indeed, these 
perturbations disappear, when relatively large time-step sizes are used (see 
Fig. 4, bottom/right). 



2 -Rls = 3 x Rs — 6 x _R g , where Rs and R g are the Schwarzschild and gravitational 
radii, respectively. 
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In addition, the semi-explicit solver has been tested in combination with 
the residual smoothing strategy. As expected, this approach accelerates the 
convergence considerably (Fig. 4: compare the plots bottom/left with the 
top/right). 

In most of the cases considered here, the time-step size is set to increase in a 
well-prescribed manner and independent of the residual. However, determining 
the size of the time step from the residual directly did not provide satisfactory 
convergence histories (Fig. 4, middle/right). 

The results obtained here indicate that the semi-explicit method is stable and 
can be applied to search for stationary solutions using large time steps, or 
cquivalcntly, CFL-numbers that are significantly larger than unity (Fig. 4, 
middle/left). 

4.2 Shock formation around black holes 

Similar to the forward facing step in CFD, a cold and dense disk has been 
placed in the innermost equatorial region: [1 < r < 10] x [—0.3 < 8 < 0.3]. We 
use the same parameters, initial and boundary condition as in the previous 
flow problem. A vanishing in- and out-flow conditions have been imposed 
at the boundaries of this disk. The gas surrounding the disk is taken to be 
inviscid, thin, hot and non-rotating. Thus, the flow configuration is similar 
to the forward facing step problem usually used for test calculations in CFD. 
The disk here serves as a barrier that forbids the gas from freely falling onto 
the black hole, and instead, it forms a curved shock front around the cold 
disk. The purpose of this test is mainly to examine the capability of the 
hierarchical scenario at employing the semi-explicit method adequately and 
enables capturing steady solution governed by strong shocks. In solving the 
HD-equations, an advection scheme of third order spatial accuracy and first 
order accurate in time has been used. The domain of calculation is sub-divided 
into 200 strongly-stretched finite volume cells in the radial direction and 60 
in horizontal direction. In Fig. 5 the configuration of the steady distributions 
of the density, temperature and the velocity field are shown. Similar to the 
calculation in the previous sub-section, the results indicate that the method 
employed is stable and converges to the sought steady solution even when a 
CFL-number of order 200 is used. However, the method converges relatively 
slowly compared to the implicit operator splitting approach, where steady 
solutions have been obtained after one thousand iterations only. 

4.3 Formation and acceleration of proton-dominated jets in active 
galaxies 

To study the mechanisms underlying jet formation around black holes, we 
have placed initially a classical accretion disk within the first 20 last stable 
radii, sandwiched by a hot and tenuous corona, and threaded by a large scale 
magnetic field. The solution procedure run as follows: 
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1. The HD-equations are solved using the IOS-approach as depicted in Stage 
I of Fig. 3. The calculations were run to cover the viscous time scale. 

2. Using the obtained results from the previous stage as starting conditions, 
Stage II of the global solution procedure is now employed to run the 
calculations for an additional viscous time scale. Here, the HD and the 
MHD equations are solved in a blockwise manner as described in Fig. 3. 
Stage III was not employed, as Alfven-waves propagation enhances the 
time-dependency of the flow even more. 

3. The final flow-configuration apparently governed by inflow and outflow 
plasmas. In general, outflows are gravitationally unbound, and there- 
fore the corresponding Bernoulli number should be positive, whereas 
negative numbers correspond to gravitationally bound flows that should 
end their motion inside the black hole. Fig. 6 shows the 2D distribu- 
tion of the Bernoulli number which obviously show the locations of the 
gravitationally-bound and unbound flows. 

4.4 The spectral energy distribution of the in- and outflow around 
the supermassive black hole of the giant elliptical galaxy M87 

The results obtained in the previous subsection are used to construct the spec- 
tral energy distribution. Therefore, the last stage of the hierarchical scenario 
is now employed in combination with Stage II. Here, the solver of Stage II is 
activated once every several dozens iterations of the RT-solver. 
In Fig. 9 we display the results of several calculations under various conditions. 
The results displayed in Fig. 7 and 8 are preliminary, as the distributions of 
the density and temperature used here are artificial, but aimed at testing the 
convergence of the RT-solver. 

5 The combined solution procedure: The hierarchical 
scenario 

In the following we describe the main steps of a possible algorithmic procedure 
for solving the combined set of MHD and the RT equations (see Stages III 
and IV of Fig. 3): 

1. Compute the RHSi and the Jacobian Ai = dLqi/dqi of each physical vari- 
able gi(= p,m,n, ...), where Lqi is the equation describing the evolution 
of variable qi. 

2. For each equation Lqi, compute the coefficient matrices Bi = dLqi/dqj, for 
which i j. This procedure applies for advection and diffusive operators 
only, though not for the source terms. 

3. Compute the coefficient matrices corresponding to the source terms only, 
i.e., Hi — dLqi/dqj, for i = 1, N and j = 1, N, and i ^ j. 
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The separation of the above-mentioned procedures is essential for enhancing 
the global efficiency of the hierarchical method. Specifically, the computation 
of each of the Bi and Hi is optional, depending on the problem in hand. For 
example, to solve the system of equations corresponding to the hydrodynam- 
ical and isothermal flow in ID efficiently, the numerical algorithm should be 
capable of calling the relevant routines only. Thus, non-relevant routines can 
be switched off almost automatically, depending on the problem in hand. In 
particular, enlarging (reducing) dimensions, incorporating additional (exclud- 
ing) variable should be algorithmically maintainable. 

Taking into account that most astrophysical flows are of multi-scale by 
nature, we think that the hierarchical solution strategy might be a promising 
approach. In the following, we describe briefly the basis of this hierarchical 
scenario applied to set of radiative MHD and the RT equations. 

1. The hierarchical approach, or equivalently the multi-stage solution proce- 
dure, is based primarily on designing the global solver in such a manner to 
achieve maximum flexibility. Specifically, the numerical algorithm should 
be capable of solving the equations sequentially, block-sequential and/or in 
a fully-coupled manner. Re-ordering and using different pre-conditioning 
should be maintainable without changing the core of the inverter. 

2. As far as vortex-free compressible, viscous and time-dependent flows are 
concerned, the implicit operator splitting approach (IOS) has been veri- 
fied to be efficient and robust. IOS is most appropriate for astrophysical 
fluid simulations, when the sought solutions depend weakly on the initial 
conditions, but strongly on the boundary conditions. The IOS-method is 
based on solving the set of equations sequentially as described in Stage I 
of Fig. 3. The convergence rate of the IOS-method may depend consider- 
ably on the order in which the equations are solved, provided the number 
of global iterations is low. 

3. The coupling between the equations can be enhanced gradually. From the 
cluster of coefficients, we may construct the Jacobians A HD and ^4 MHD ; 
which correspond to the set of HD and MHD equations (see Stage II/Fig. 
3). Algorithmically, this procedure is basically a sort of re-ordering and re- 
organizing of the coefficients, and does not require an extensive program- 
ming. As in the previous step, the order in which the equations are solved 
may affect both its convergence rate and efficiency. Here, a special care 
should be given to assure that the inclusion of coefficients corresponding 
to the source terms does not enlarge the band width of A HD and ^4 MHD . 
Test calculations have shown that careful ordering of the HD-equations 
may re duce the computational costs devoted for matrix inversion by 75% 
i HU.Ilh . Furthermore, it has been verified that several equations can still 
be separated and solved sequentially. Namely, the Possion equation for 
modelling self-gravity as well as the angular momentum equation accept 
partial decoupling from the rest of equations, provided the flow is axi- 
symmetric. 
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4. Using the solutions obtained in stage II as initial conditions, we may solve 
the whole set of HD and MHD equations as a single set of coupled equa- 
tions. The resulting Jacobian is highly sparse, for which pre-conditioned 
Krylov sub-iterative methods are highly appropriate. 

5. By iterating over Stage II and IV, we can be sure that the resulting solu- 
tion is reasonably close to sought quasi-stationary or steady solutions for 
the radiative MHD and radiative transfer equations. This is a consequence 
of: 

a) The radiative intensity in the high density regions, where the optical 
thickness is large, is isotropic and coincides with black-body emission. 
Therefore, the intensity obtained by solving the zero moment of the 
radiation field is sufficiently accurate in this regime. 

b) The radiative intensity obtained by solving the RT-equation in opti- 
cally thin regions may differ considerably from that obtained using 
the gray approximation. However, radiation in such regions have neg- 
ligible power and they may hardly affect the dynamics of the flow. 
Consequently, the following solution method may be proposed: 

• The numerical values of the variables obtained in Stage III are 
used as initial conditions for calculating the non-gray and time- 
dependent radiative intensity. 

• The mean- value of the frequency-dependent intensity is computed 
and subsequently used as initial condition for the radiative MHD 
equations. 

• To avoid extensive computational costs, it is suggested to solve for 
I v every 10, or 20 time-steps. However, since the radiative time- 
scale is extremely short compared to the hydrodynamical time 
scale, it is much more reasonable to solve for the time-independent 
intensity. 

6 Summary 

In this paper we have presented the hierarchical scenario for solving the set of 
radiative MHD equations and the 5D axi-symmetric radiative transfer equa- 
tion. 

The main features of this scenario are as follows: 

1. The global efficiency can be enhanced, depending on the optimal architec- 
ture of the global solver. Specifically, the algorithmic structure should be 
sufficiently flexible, so that scalar or set of equations in arbitrary dimen- 
sions, different accuracies and using the appropriate pre-conditionings can 
be solved with a reasonable efficiency. 

2. Robustness is monitored through employing a variety of solution pro- 
cedures. Depending on the particular features of the problem considered, 
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several stages of implicitness may be used, depending on the number of co- 
efficients used for constructing the coefficient matrix. In particular, start- 
ing with a purely explicit time-stepping scheme, the algorithm should be 
capable of modifying the scheme into a fully implicit method dynamically. 

3. For implicit calculations, the hierarchical algorithm relies on using a vari- 
ety of preconditioning for accelerating convergence. For example, for mod- 
elling weakly incompressible flows, it has been verified that the "Approx- 
imate Factorization Method" as pre-conditioning yields a larger conver- 
gence rate than the "Alternating Directional Implicit" or the "Line Gauss- 
Seidel" methods. However, the latter preconditionings provide faster con- 
vergence if the flow is compressible and advection-dominatcd. Therefore, 
depending on the problem in hand, the algorithm should be capable of 
employing the appropriate preconditioning at least in an explicit-adaptive 
manner. 

4. The hierarchical algorithm is capable of solving the angle-averaged time- 
dependent radiation transfer equation, taking into account the Kompa- 
neets operator for modelling up-scattering of soft photons by hot electrons 
in magnetized plasmas. 

We note, however, that the assumption of isotropic radiative intensity may 
break down if the flow is relativistic and contains regions of significantly 
different optical depths. Therefore, in the near future we intend to modify 
the RT-solver to enable modelling the motions of ultra-relativistic plasmas 
in the vicinities Kerr and Schwarzschild black holes. 

5. The algorithm includes a procedure that allows solving the zero-moment 
MHD equations partially/loosely coupled with the radiation transfer equa- 
tions. The latter coupling can be significantly enhanced through paral- 
lelization on powerful machines. 

Finally, we have shown that the hierarchical algorithm presented here can 
be applied to study the mechanisms underlying the formation, launching and 
acceleration of jets in AGNs and quasars, though serious numerical and phys- 
ical modifications are still required. 

References 

[BA91] Balbus, S., Hawley, J., 1991, "A powerful local shear instability 

in weakly magnetized disks. I - Linear analysis. II - Nonlinear 

evolution", ApJ, 376, 214 
[BW78] Beam, R.M., Warming, R.F., 1978, "An implicit factorized scheme 

for compressible Navier-Stokes equations" ,AIAA, 16, 393 
[BR01] Brandt, A., 2001, "Textbook Multi-grids" , in Multigrid, ed.: Trot- 

tenberg, U., Oosterlee, C., Schuller, A., Acad. Press, London 
[D098] Dongarra, I.S., Duff, D.C., Sorcnscn, H.A., van dcr Vorst, 1998, 

"Num. Linear Alg. for High-Pcrformancc Computers", SIAM J. 

Scicnt. Comput., 20, 94 



26 Hujeirat, A. 



[FE72] Felten, J.E., & Rees, M.J., 1972, "Transfer effects on lines and 
continuum in optically thick sources", A&A, 21, 139 

[FL88] Fletcher, C.A.J. , 1988, 'Computational Techniques for Fluid Dy- 
namics', Vol, I and II, Springer- Verlag 

[FO00] Font, J. A. 2000, "Numerical Hydrodynamics in General Relativ- 
ity", Living Rev. Relativity, 3, 2 

[FR00] Fryxell, B., Olson, K., Ricker, P., et al., 2000, "FLASH: An Adap- 
tive Mesh Hydrodynamics Code for Modeling Astrophysical Ther- 
monuclear Flashes", ApJS, 131, 273 

[GA03] Gammic, C.F., McKinney, J.C., Toth, G., 2003, "HARM: A Nu- 
merical Scheme for General Relativistic Magnctohydrodynamics" , 
ApJ, 589, 444 

[HI90] Hirsch, C, 1990, 'Num. Computation of Internal and External 
Flows', Vol, I, and II, John Wiley & Sons, New York 

[HUJ0] Hujeirat, A., Papaloizou, J. CP., 1998, "Shock formation in ac- 
cretion columns - a 2D radiative MHD approach" ,A&A, 340, 593 

[HUJ1] Hujeirat, A., Rannacher,R., 2001, "On the efficiency and robust- 
ness of implicit methods in computational astrophysics" , NewAR, 
45, 425 

[HUJ2] Hujeirat, A., Camenzind, M., Livio, M., 2002, "Ion-dominated 
plasma and the origin of jets in quasars" , A&A, 394, L9 

[HUJ3] Hujeirat, A., Camenzind, M., Burkert, A., 2002b, "Comptoniza- 
tion and synchrotron emission in 2D accretion flows. I. A new 
numerical solver for the Kompaneets equation", A&A, 386 , 757 

[HUJ5] Hujeirat, A., Livio, M., Camenzind, M., Burkert, A., 2003, "A 
model for the jet-disk connection in BH accreting systems" , A&A, 
408, 415 

[HUJ6] Hujeirat, A., Blandford, R.D., 2004, "A model for electromag- 
netic extraction of rotational energy and formation of accretion- 
powered jets in radio galaxies", A&A, 416, 423 

[HUJ8] Hujeirat, A., 2004, "A method for enhancing the stability and 
robustness of explicit schemes in CFD" , in press, New Astronomy 
Reviews. 

[KA76] Katz, J. A., 1976, "Nonrelativistic Compton scattering and models 
of quasars", ApJ, 206, 910 
[IL72] Iilarinov, A.F., & Sunyaev, R.A., 1972, "Compton scattering by 
thermal electrons in X-ray sources "Soviet Astr. -AJ, 16, 45 

[KL89] Kley, W., 1989, "Radiation hydrodynamics of the boundary layer 
in accretion disks. I - Numerical methods", A&A, 208, 98 

[K099] Koide, S., Shibata, K., & Kudoh, T. 1999, "Relativistic Jet For- 
mation from Black Hole Magnetized Accretion Disks: Method, 
Tests, and Applications of a General Relativistic Magnetohydro- 
dynamic Numerical Code", ApJ, 522, 727 



Hierarchical scenario: Radiative MHD solvers 



27 



[KO02] Koide, S., Shibata, K., Kudoh, T., & Meier, D. L. 2002, "Extrac- 
tion of Black Hole Rotational Energy by a Magnetic Field and 
the Formation of Relativistic Jets", Science, 195, 1688 

[K099] Komissarov, S. S. 1999, "A Godunov-type scheme for relativistic 

magnetohydrodynamics" , MNRAS, 303, 343 
[LE81] Levermore, CD., & Pomraning, G.C., 1981, "A flux- limited dif- 
fusion theory", ApJ, 248, 321 

[MA96] Mahadevan, R., & Narayan, R., Yi, I., 1996, " Harmony of elec- 
trons: Cyclotron and Synchrotron emission by thermal electrons 
in magnetic fields", ApJ, 465, 327 

[MA85] MacCormack, R.W., 1985, "Current status of numerical solutions 
of Navicr-Stokes equations", AIAA, Paper 81-0110 

[MA99] Marti', J.M., Muller, E., 1999, "Numerical hydrodynamics in spe- 
cial relativity", Living Rev. Relativity, 2, 3 

[ME01] Meier, D.L., Koide, S., & Uchida, Y. 2001, "Magnetohydrody- 
namic Production of Relativistic Jets", Science, 291, 84 

[ME03] Meier, D., 2003, "The theory and simulation of relativistic jet 
formation: towards a unified model for micro- and macroquasars" , 
NewAR, 47, 667 

[MI86] Mihalas, D., Mihalas, B.W., 1984, "Foundations of radiation hy- 
drodynamics" , Oxford University Press, NY, (MM) 

[OU97] Ouyed, R., Pudritz, R., 1997, "Numerical simulation of astrophys- 
ical jets from Kcplcrian disks. II. episodic outflows", ApJ, 484, 
794 

[PA80] Payne, D.G., 1980, "Time-dependent Comptonization - X-ray re- 
verberations" , ApJ, 237, 951 

[RY79] Rybiki, G.B., & Lightman, A. P., 1979, Radiation processes, 
Wiley-Interscience Publication 

[SA00] Saad, Y., van der Vorst, 2000, "Iterative solution of linear systems 
in the 20-th century", J. of Comp. and Appl. Math., 123, 1 

[SH76] Shapiro, S.L., Lightman A.P., k Eardley, D.M., "A two- 
temperature accretion disk model for Cygnus X-l structure and 
spectrum", 1976, ApJ, 204, 187 

[ST92] Stone, J.M., & Norman, M., 1992, "ZEUS-2D: A radiation mag- 
netohydrodynamics code for astrophysical flows in two space di- 
mensions. I - The hydrodynamic algorithms and tests.", ApJS, 
80, 791 

[T098] Toth, Keppens, R., Botchev, M.A., 1998, "Implicit and semi- 
implicit schemes in the Versatile Advection Code: numerical 
tests", A&A, 332, 1159 

[TR01] Trottenberg, U., 2001, in Multigrid, ed.: Trottcnbcrg, U., Ooster- 
lee, C, Schuller, A., Acad. Press, London 

[UC99] Uchida, Y., Nakamura, M., Hirose, S., Uemura, S., "Magnetody- 
namic formation of jets in accretion process of magnetized mass 
onto the central gravitator", Ap&SS, 264, 195 



Hujeirat, A. 



[VI03] De Villiers, J.-R, & Hawley, J.F., 2003, "A Numerical Method for 
General Relativistic Magnetohydrodynamics" , ApJ, 589, 458 

[ZI98] Ziegler, U., 1998, "NIRVANA+: An adaptive mesh refinement 
code for gas dynamics and MHD", Comp. Phys. Comm., 109, 
142 



This figure "Fig5.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/math-ph/04 1 005 5 v 1 



This figure "Fig6.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/math-ph/04 1 005 5 v 1 



This figure "Fig9a.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/math-ph/04 1 005 5 v 1 



